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A TUTORIAL ON INVERSE PROBLEMS FOR ANOMALOUS DIFFUSION PROCESSES 
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Abstract. Over the last two decades, anomalous diffusion processes in which the mean squares variance grows 
slower or faster than that in a Gaussian process have found many applications. At a macroscopic level, these 
processes are adequately described by fractional differential equations, which involves fractional derivatives in 
time or/and space. The fractional derivatives describe either history mechanism or long range interactions 
of particle motions at a microscopic level. The new physics can change dramatically the behavior of the 
forward problems. For example, the solution operator of the time fractional diffusion diffusion equation has 
only limited smoothing property, whereas the solution for the space fractional diffusion equation may contain 
weakly singularity. Naturally one expects that the new physics will impact related inverse problems in terms 
of uniqueness, stability, and degree of ill-posedness. The last aspect is especially important from a practical 
point of view, i.e., stably reconstructing the quantities of interest. 

In this paper, we employ a formal analytic and numerical way, especially the two-parameter Mittag-Leffler 
function and singular value decomposition, to examine the degree of ill-posedness of several “classical” inverse 
problems for fractional differential equations involving a Djrbashian-Caputo fractional derivative in either time 
or space, which represent the fractional analogues of that for classical integral order differential equations. We 
discuss four inverse problems, i.e., backward fractional diffusion, sideways problem, inverse source problem and 
inverse potential problem for time fractional diffusion, and inverse Sturm-Liouville problem, Cauchy problem, 
backward fractional diffusion and sideways problem for space fractional diffusion. It is found that contrary 
to the wide belief, the influence of anomalous diffusion on the degree of ill-posedness is not definitive: it can 
either significantly improve or worsen the conditioning of related inverse problems, depending crucially on 
the specific type of given data and quantity of interest. Further, the study exhibits distinct new features of 
“fractional” inverse problems, and a partial list of surprising observations is given below. 

(a) Classical backward diffusion is exponentially ill-posed, whereas time fractional backward diffusion is only 
mildly ill-posed in the sense of norms on the domain and range spaces. However, this does not imply 
that the latter always allows a more effective reconstruction. 

(b) Theoretically, the time fractional sideways problem is severely ill-posed like its classical counterpart, but 
numerically can be nearly well-posed. 

(c) The classical Sturm-Liouville problem requires two pieces of spectral data to uniquely determine a general 
potential, but in the fractional case, one single Dirichlet spectrum may suffice. 

(d) The space fractional sideways problem can be far more or far less ill-posed than the classical counterpart, 
depending on the location of the lateral Cauchy data. 

In many cases, the precise mechanism of these surprising observations is unclear, and awaits further analyti¬ 
cal and numerical exploration, which requires new mathematical tools and ingenuities. Further, our findings 
indicate fractional diffusion inverse problems also provide an excellent case study in the differences between 
theoretical ill-conditioning involving domain and range norms and the numerical analysis of a finite-dimensional 
reconstruction procedure. Throughout we will also describe known analytical and numerical results in the lit¬ 
erature. 

Keywords: fractional inverse problem; fractional differential equation; anomalous diffusion; Djrbashian- 
Caputo fractional derivative; Mittag-Leffler function. 


1. Introduction 

Diffusion is one of the most prominent transport mechanisms found in nature. At a microscopic level, it is 
related to the random motion of individual particles, and the use of the Laplace operator and the first-order 
derivative in the canonical diffusion model rests on a Gaussian process assumption on the particle motion, after 
Albert Einstein’s groundbreaking work [23]. Over the last two decades a large body of literature has shown that 
anomalous diffusion models in which the mean square variance grows faster (superdiffusion) or slower (subdiffu¬ 
sion) than that in a Gaussian process under certain circumstances can offer a superior fit to experimental data 
(see the comprehensive reviews ED ED] □ tmi for physical background and practical applications). For example, 
anomalous diffusion is often observed in materials with memory, e.g., viscoelastic materials, and heterogeneous 
media, such as soil, heterogeneous aquifer, and underground fluid flow. At a microscopic level, the subdiffusion 
process can be described by a continuous time random walk [75], where the waiting time of particle jumps follows 
some heavy tailed distribution, whereas the superdiffusion process can be described by Levy flights or Levy walk, 
where the length of particle jumps follows some heavy tailed distribution, reflecting the long-range interactions 
among particles. Following the aforementioned micro-macro correspondence, the macroscopic counterpart of a 
continuous time random walk is a differential equation with a fractional derivative in time, and that for a Levy 
flight is a differential equation with a fractional derivative in space. We will refer to these two cases as time 
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fractional diffusion and space fractional diffusion, respectively, and it is generically called a fractional derivative 
equation (FDE) below. In general the fractional derivative can appear in both time and space variables. 

Next we give the mathematical model. Let f 1 = (0,1) be the unit interval. Then a general one-dimensional 
FDE is given by 

(1.1) dtU - qD%u + qu = / (x,t) G ft x (0,T), 

where T > 0 is a fixed time, and it is equipped with suitable boundary and initial conditions. The fractional 
orders a £ (0,1) and /3 £ (1, 2) are related to the parameters specifying the large-time behavior of the waiting-time 
distribution or long-range behavior of the particle jump distribution. For example, in hydrological studies, the 
parameter /3 is used to characterize the heterogeneity of porous medium m In theory, these parameters can be 
determined from the underlying stochastic model, but often in practice, they are determined from experimental 
data [34l [351 [ 6 T] . The notation d “ = is the Djrbashian-Caputo derivative operator of order a £ (0,1) in the 
time variable t, and qD£ denotes the Djrbashian-Caputo derivative of order /? £ (1,2) in the space variable x. 
For a real number n — l<j<n, n£N, and / £ H n ( 0,1), the left-sided Djrbashian-Caputo derivative oD2 f of 
order 7 is defined by G20 pp- 9i] 

(1-2) [lD2f = r| , 1 . f\x - s) n - 1 -y^ ) (s)ds, 

r(«-7) Jo 

where T(z) denotes Euler’s Gamma function defined by 

POO 

r(z) = / s z ~ 1 e~ s ds, K(z) > 0 . 

Jo 

The Djrbashian-Caputo derivative was first introduced by Armenian mathematician Mkhitar M. Djrbashian for 
studies on space of analytical functions and integral transforms in 1960s (see [2ll l20l [19] for surveys on related 
works). Italian geophysicist Michele Caputo independently proposed the use of the derivative for modeling the 
dynamics of viscoelastic materials in 1967 [TO]. We note that there are several alternative (and different) definitions 
of fractional derivatives, notably the Riemann-Liouville fractional derivative, which formally is obtained from 
ED by interchanging the order of integration and differentiation, i.e., the left-sided Riemann-Liouville fractional 
derivative / of order 7 £ (n — 1, n), n £ N, is defined by [531 pp. 70] 

jn 1 rx 

Rplfm-?——-—- (x-s) n ~ 1 -' y f(s)ds. 

0 J dx n r(n — 7 ) J 0 v J w 

In this work, we shall focus mostly on the Djrbashian-Caputo derivative since it allows a convenient treatment of 
the boundary and initial conditions. 

Under certain regularity assumption on the functions, with an integer order 7 , the Djrbashian-Caputo and 
Riemann-Liouville derivatives both recover the usual integral order derivative (see for example 1791 pp. 100] for the 
Djrbashian-Caputo case). For example, with a = 1 and /3 = 2, the Djrbashian-Caputo fractional derivatives <9fu 
and o'D" u coincide with the usual first- and second-order derivatives and , respectively, for which the model 
El recovers the standard one-dimensional diffusion equation, and thus generally the model El is regarded as 
a fractional counterpart. The Djrbashian-Caputo derivative (and many others) is an integro-differential operator, 
and thus it is nonlocal in nature. As a consequence, many useful rules, e.g., product rule and integration by parts, 
from PDEs are either invalid or require significant modifications. The nonlocality underlies most analytical and 
numerical challenges associated with the model ED- It significantly complicates the mathematical and numerical 
analysis of the model, including relevant inverse problems. 

In a fractional model, there are a number of parameters, e.g., fractional order(s), diffusion and potential coef¬ 
ficients (when using a second-order elliptic operator in space), initial condition, source term, boundary conditions 
and domain geometry, that cannot be measured/specified directly, and have to be inferred indirectly from mea¬ 
sured data. Typically, the data is the forward solution restricted to either the boundary or the interior of the 
physical domain. This gives rise to a large variety of inverse problems for FDEs, which have started to attract 
much attention in recent years, since the pioneering work M- An interesting question is how the nonlocal physics 
behind anomalous diffusion processes will influence the behavior of related inverse problems, e.g., uniqueness, sta¬ 
bility, and the degree of ill-posedness. The degree of ill-posedness is especially important for developing practical 
numerical reconstruction procedures. There is a now well known example of backward fractional diffusion, i.e., 
recovering the initial condition in a time fractional diffusion equation from the final time data, which is only mildly 
ill-posed, instead of severely ill-posed for the classical backward diffusion problem. I 11 some sense, this example 
has led to the belief that “fractionalizing” inverse problems can always mitigate the degree of ill-posedness, and 
thus allows a better chance of an accurate numerical reconstruction. 
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In this paper, we examine the degree of ill-posedness of “fractional” inverse problems from a formal analytic and 
numerical point of view, and contrast their numerical stability properties with their classical, that is, the Gaussian 
diffusion counterparts, for which there are many deep analytical results |84l 1401 1301 . Specifically, we revisit 
a number of “classical” inverse problems for the FDEs, e.g., the backward diffusion problem, sideways diffusion 
problem and inverse source problem, and numerically exhibit their degree of ill-posedness. These examples indicate 
that the answer to the aforementioned question is not definitive: it depends crucially on the type (unknown and 
data) of the inverse problem we look at, and the nonlocality of the problem (fractional derivative) can either 
greatly improve or worsen the degree of ill-posedness. 

The mathematical theory of inverse problems for FDEs is still in its infancy, and thus in this work, we only 
discuss the topic formally to give a flavor of inverse problems for FDEs - our goal is to give insight rather than 
to pursue an in-depth analysis. The technical developments that are available we leave to the references cited. In 
addition, known theoretical results and computational techniques in the literature will be briefly described, which 
however are not meant to be exhaustive. The rest of the paper is organized as follows. In Section[2]we review two 
special functions, i.e., Mittag-Leffler function and Wright function, and their basic properties. The Mittag-LefHer 
function plays an extremely important role in understanding anomalous diffusion processes. We also recall the 
basic tool - singular value decomposition - for analyzing discrete inverse problems. Then in Section [3] we study 
several inverse problems for FDEs with a time fractional derivative, including backward diffusion, inverse source 
problem, sideways problem and inverse potential problem. In Section[I]we consider inverse problems for FDEs with 
a space fractional derivative, including the inverse Sturm-Liouville problem, Cauchy problem, backward diffusion 
and sideways problem. In the appendices, we give the implementation details of the computational methods 
for solving the time- and space fractional differential equations. These methods are employed throughout for 
computing the forward map (unknown-to-measurement map) so as to gain insight into related inverse problems. 
Throughout the notation c, with or without a subscript, denote a generic constant, which may differ at different 
occurrences, but it is always independent of the unknown of interest. 


2. Preliminaries 

We recall two important special functions, Mittag-Leffler function and Wright function, and one useful tool for 
analyzing discrete ill-posed problems, singular value decomposition. 


2.1. Mittag-Leffler function. We shall use extensively the two-parameter Mittag-Leffler function E a ^(z) (with 
a > 0 and /3 £ R) defined by 1531 equation (1.8.17), pp. 40] 


( 2 . 1 ) 


-£'q:, ( 5 (-Z) 


oo 


E 


r (ka + /3) 


z £ €. 


This function with /3 = 1 was first introduced by Gosta Mittag-Leffler in 1903 m and then generalized by others 
UDE. It can be verified directly that 


Ei,i(z) = e z , E 2 ,i(z) = cosh yfz, E 2 ,i(z) 


sinh y/z 


Hence it represents a generalization of the exponential function in that E\,\(z) = e z . The Mittag-Leffler function 
E a ,p(z) is an entire function of z with order a” 1 and type 1 [52, pp. 40]. Further, the function E a ,i{—t) is 
completely monotone on the positive real axis R + 1821 . and thus it is positive on R + ; see also [90] for extension 
to the two-parameter Mittag-Leffler function E a ,p(z). It appears in the solution representation for FDEs: The 
functions E a p(—Xt a ) and t a ~ 1 E a , a (—\t a ) appear in the kernel of the time fractional diffusion problem with 
initial data and the right hand side, respectively, and also are eigenfunctions to the fractional Sturm-Liouville 
problem with a zero potential, cf. Section [4M~1 

In our discussions, the asymptotic behavior of the function E a ^(z) will play a crucial role. It satisfies the fol¬ 
lowing exponential asymptotics 1531 pp. 43, equations (2.8.17) and (2.8.18)], which was first derived by Djrbashian 
121] , and refined by many researchers 110511801 . 


Lemma 2.1. Let a £ (0,2), f3 £ R, and fi £ (an/2, min(7r, an)). 

1 N 
1 n-fi 


E a A*) = -z (1 - WA V 1/ “ - E '4 + 0(- 

a J F(p — ak) z“ a 


E a Az) = - E r(/g — ak) ~zX + ° ( ^TT ) wUh |z| 


Then for N £ N 

IT+i) with \z\ -> oo, |arg(z)| < n, 
oo , At < |arg(«)| < n. 
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From these asymptotics, the Mittag-Leffler function E a ^(z) decays only linearly on the negative real axis 
R“, which is much slower than the exponential decay for the exponential function e z . However, on the positive 
real axis R + , it grows exponentially, and the growth rate increases with the fractional order 0 < a < 2. To 
illustrate the distinct feature, we plot the functions E a ,i{—Tv 2 t a ) and t a ~ 1 E c , }a (— 7r 2 t Q ) in Fig. [T] for several 
different a values, where A = 7r 2 is the first Dirichlet eigenvalue of the negative Laplacian on the unit interval 
fl = (0,1); see Appendix I A. II for further details on the computation of the Mittag-Leffler function. Fig. 1(a) 
can be viewed as the time evolution of u(l/2,t), where <3“u — u xx = 0 with u(0, t) = u(l,f) = 0, and initial 
data uo{x) = sin7ra; (the lowest Fourier eigenmode). The slow decay behavior at large time is clearly observed. 
In particular, at t = 1, the function -E q ,i(— iv 2 t) still takes values distinctly away from zero for any 0 < a < 1, 
whereas the exponential function e almost vanishes identically. In contrast, for t close to zero, the picture is 
reversed: the Mittag-Leffler function E a ^{—n t) decays much faster than the exponential function e -7r . The 
drastically different behavior of the function f? a ,i(— z), in comparison with the exponential function e _z , explains 
many unusual phenomena with inverse problems for FDEs to be described below. According to the exponential 
asymptotics, the function E a , a (z) decays faster on the negative real axis R“, since l/r(0) = 0, i.e., the first term 
in the expansion vanishes. This is confirmed numerically in Fig. Eb). Even though not shown, it is noted that 
the function E a:a (z) decays only quadratically on the negative real axis R“ for a € (0,1) or a £ (1, 2), which is 
asymptotically much slower than the exponential decay. 




(a) E aA {-n 2 t a ) (b) r- 1 E Q , a (- ttV) 

Figure 1. The profiles of Mittag-Leffler functions (a) E a: i(—n 2 t a ) and (b) t a ~ 1 E a ^ a {—n 2 t a ). 
The value 7r 2 is the first Dirichlet eigenvalue of the negative Laplacian on the interval (0,1). 


The distribution of zeros of the Mittag-Leffer function E a ,p(z) is of immense interest, especially in the related 
Sturm-Liouville problem; see Section I4T1 below. The case of /3 = 1 was first studied by Wiman [104] . It was 
revisited by Djrbashian m, and many deep results were derived, especially for the case of a = 2. There are 
many further refinements [93]; see [83] for an updated account. 

2.2. Wright function. The Wright function W P:IJ ,(z) is defined by 110611 1071 

OO £ 

WpAz) = ^] m Z pk+»y /x ’ p6K ’ p >- 1 ’ * GC - 

This is an entire function of order 1/(1 + p) [301 Theorem 2.4.1], It was first introduced in connection with a 
problem in number theory by Edward M. Wright, and revived in recent years since it appears as the fundamental 
solution for FDEs [69] . The Wright function W p ^{z) has the following the asymptotic expansion in one sector 
containing the negative real axis r- Theorem 3.2]. Like before, the exponential asymptotics can be used to 
deduce the distribution of its zeros [66| . 
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Lemma 2 . 2 . Let — 1 < p < 0, y = —z, arg (z) < tv, —tv < arg (y) < tv, |arg(y)| < min(37r(l + p)/2, tv) — e, e > 0. 
Then 


W^z) = Y 


- V 1/2~H 0 -Y 


Y A m Y- m + 0(Y- 


Y —> oo, 


where Y = (1 + p)((—p) P yfi^ 1+ A and the coefficients A m , m = 0, 1,.. . are defined by the asymptotic expansion 


T(1 -n-pt) 


2tt(— p)-^(l + p )(i+/>)(t+i)r(t + 1) 


•£ 


(-i ) m A„ 


= 0 r((i + p)t + p + \ + m ) 
+ o‘ 1 


v r((i + p)£-f-/3 + ^ + m) y 

valid for arg(t), arg(— pt), and arg(l — p — pt) all lying between —n and n and t tending to infinity. 


The Wright function W P}IJ ,(z), — 1 < p < 0, decays exponentially on the negative real axis R - , in a manner 
similar to the exponential function e z , but at a different decay rate. Its special role in fractional calculus is 
underscored by the fact that it forms the free-space fundamental solution K a (x,t) to the one-dimensional time 
fractional diffusion equation [69] by 

(2.2) K a (x,t) = W_ fil _ f (-M/r /2 ) . 


The multidimensional case is more complex and involves further special functions, in particular, the Fox H function 
[9111541 . For a = 1, the formula (12.211 recovers the familiar free-space fundamental solution for the one-dimensional 
heat equation, i.e., 


K(x, t) = 


1 


2\fivi 

which is a Gaussian distribution in x for any t > 0. In the fractional case, the fundamental solution K a (x,t) 
exhibits quite different behavior than the heat kernel. To see this, we show the profile of K a (x,t) in Fig. [2] for 
several a values; see Appendix ro for a brief description of the computational details. For any 0 < a < 1, 
K a (x,t) decays slower at a polynomial rate as the argument tends to infinity, i.e., having a long tail, 

when compared with the Gaussian density. The long tail profile is one of distinct features of slow diffusion [5]. 
Further, for any a < 1, the profile is only continuous but not differentiable at x = 0. The kink at the origin 
implies that the solution operator to time fractional diffusion may only have a limited smoothing property. 




X X 

(a) t = 0.1 (b) t = 1 


Figure 2. The profile of the fundamental solution K a (x,t) at (a) t = 0.1 and (b) t = 1. 
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2.3. Singular value decomposition. We shall follow the well-established practice in the inverse problem com¬ 
munity, i.e., using the singular value decomposition, as the main tool for numerically analyzing the problem 
behavior [32] , Specifically, we shall numerically compute the forward map F, and analyze its behavior to gain 
insights into the inverse problem. Given a matrix A £ R rixm , its singular value decomposition is given by 

A = UEV\ 

where U = [ui ... u n ] £ R™ X71 and V = [vi ... v m ] £ R mxm are column orthonormal matrices, and E £ 
R™ xm _ c jiag( (Jl! .... CTr , ) 0,..., 0) is a diagonal matrix, with the diagonal elements {cr^} being nonnegative and 
listed in a descending order or > ... > ay > 0, and r being the (numerical) rank of the matrix A. The diagonal 
element a; is known as the ith singular value, and the corresponding columns of U and V, i.e., m and Vi, are 
called the left and right singular vectors, respectively. 

One simple measure of the conditioning of a linear inverse problem Ax = b is the condition number cond(A), 
which is defined as the ratio of the largest to the smallest nonzero singular value, i.e., 

cond(A) = or/ov. 

In particular, if the condition number is small, then the data error will not be amplified much. In the case of a 
large condition number, the issue is more delicate: it may or may not amplify the data perturbation greatly. A 
more complete picture is provided by the singular value spectrum (or, 02 ,.. . , ay). Especially, a singular value 
spectrum gradually decaying to zero without a clear gap is characteristic of many discrete ill-posed problems, 
which is reminiscent of the spectral behavior of compact operators. We shall adopt these simple tools to analyze 
related inverse problems below. 

In addition, using singular value decomposition and regularization techniques, e.g. Tikhonov regularization or 
truncated singular value decomposition, one can conveniently obtain numerical reconstructions, even though this 
might not be the most efficient way to do so. However, we shall not delve into the extremely important question of 
practical reconstructions, since it relies heavily on a priori knowledge on the sought for solution and the statistical 
nature (Gaussian, Poisson, Laplace ...) of the contaminating noise in the data, which will depend very much on 
the specific application. We refer interested readers to the monographs [2611921141] and the survey [49] for updated 
accounts on regularization methods for constructing stable reconstructing procedures and efficient computational 
techniques. We will also briefly mention below related works on the application of regularization techniques to 
inverse problems for FDEs. 


3. Inverse problems for time fractional diffusion 

In this section, we consider several model inverse problems for the following one-dimensional time fractional 
diffusion equation on the unit interval fl = (0,1): 

(3.1) dtu—u xx + qu = f infix (0,T], 

with the fractional order a £ (0, 1), the initial condition u(0) = v and suitable boundary conditions. Although 
we consider only the one-dimensional model, the analysis and computation in this part can be extended into 
the general multi-dimensional case, upon suitable modifications. Recall that dfu denotes the Djrbashian-Caputo 
fractional derivative of order a with respect to time t. For a = 1, the fractional derivative coincides with the 
usual first-order derivative u' , and accordingly, the model EH reduces to the classical diffusion equation. Hence 
it is natural to compare inverse problems for the model eh with that for the standard diffusion equation. We 
shall discuss the following four inverse problems, i.e., the backward problem, sideways problem, inverse source 
problem and inverse potential problem. In the first three cases, we shall assume a zero potential q = 0. We will 
aslo discuss the solution of an inverse coefficient problem using fractional calculus. 

3.1. Backward fractional diffusion. First we consider the time fractional backward diffusion. By the linearity 
of the inverse problem, we may assume that equation EH) is prescribed with a homogeneous Dirichlet boundary 
condition, i.e., u = 0 at x — 0, 1, and the initial condition u(0) = v. Then the inverse problem reads: given the 
final time data g = u(T), find the initial condition v. It arises in, for example, the determination of a stationary 
contaminant source in underground fluid flow. 

To gain insight, we apply the separation of variables. Let {(Ay, (pj)}, with A j = (jV) 2 and 4>j = %/2sin jnx, be 
the Dirichlet eigenpairs of the negative Laplacian on the interval fi. The eigenfunctions {4>j} form an orthonormal 
basis of the L 2 (fi) space. Then using the Mittag-Leffler function E a ^(z) defined in (12. Ill , the solution u to equation 
EH can be expressed as 

OO 

u(x,t) = ^E aA {-\ j t a ){y,(l> j )(l> j {x). 

0 = 1 
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Therefore, the final time data g = u(T) is given by 


3 = 1 

It follows directly that the initial data v is formally given by 


= E 




rt E a} 1 (—XjT a ) 


Since the function E a ^{—t) is completely monotone on the positive real axis R + [82] for any a £ (0,1], the 
denominator in the representation does not vanish. In case of a = 1, the formula reduces to the familiar 
expression 

oo 

v = ( 3 , 

i=i 

This formula shows clearly the well-known severely ill-posed nature of the backward diffusion problem: the 
perturbation in the jth Fourier mode (g,<f>j) of the (noisy) data g is amplified by an exponentially growing factor 
e XjT , which can be astronomically large, even for a very small index j, if the terminal time T is not very small. 
Hence it is always severely ill-conditioned and we must multiply the jth Fourier mode of the data g by a factor 
e XjT in order to recover the corresponding mode of the initial data v 

In the fractional case, by Lemma |2.II the Mittag-Leffler function E a: i(z) decays only linearly on the negative 
real axis R _ , and thus the multiplier 1/E a ,i(— XjT a ) grows only linearly in Xj, i.e., 1/E a ,i(— XjT a ) ~ A j, which is 
very mild compared to the exponential growth e AjT for the case a = 1, and thus the fractional case is only mildly 
ill-posed. Roughly, the jth Fourier mode of the initial data v now equals the jth mode of the data g multiplied 
by Xj. More precisely, it amounts to the loss of two spatial derivatives |88l Theorem 4.1] 

(3-2) IMIr, 3 (n) < c ll u Cnilff 2 (n)- 

Intuitively, the history mechanism of the anomalous diffusion process retains the complete dynamics of the physical 
process, including the initial data, and thus it is much easier to go backwards to the initial state v. This is in 
sharp contrast to the classical diffusion, which has only a short memory and loses track of the preceding states 
quickly. This result has become quite well-known in the inverse problems community and has contributed to a 
belief that “inverse problems for FDEs are less ill-conditioned than their classical counterparts” - throughout this 
paper we will see that this conclusion as a general statement can be quite far from the truth. 

Does this mean that for all terminal time T the fractional case is always less ill-posed than the classical one? 
The answer is yes, in the sense of the norm on the data space in which the data g lies. Does this mean that from 
a computational stability standpoint that one can always solve the backward fractional problem more effectively 
than for the classical case? The answer is no, and the difference can be substantial. To illustrate the point, let J 
be the highest frequency mode required of the initial data v and assume that we believe we are able to multiply 
the first J modes gj = (g, 4>j ), j — 1, 2,..., J, by a factor no larger than M (which roughly assumes that the noise 
levels in both cases are comparable). By the monotonicity of the function E at i(—t) in t, it suffices to examine the 
Jth mode. For the heat equation vj := (v,<j>j) = e XjT gj and provided that T = Tj < Aj/logM this is feasible. 
For a fixed J, if T£ denotes the point where 

e~ XjT °‘ = E aA (-XjT*), 

then in the fractional case for T < T* the growth factor on gj will exceed M for any T <T*. In Tabic Q] below, 
we present the critical value T* for several values of the fractional order a and the maximum number of modes 
J. The numbers in the table are very telling. For example, for the case J = 5, a = 1/4 and T = 0.02 (which is 
approximately one half the value of T*), the growth factor is about 1.6 for the heat equation but about 113 for 
the fractional case. With J = 10 and a = 1/4 and T = T* the growth factor is around 336. If T = T*/10 then 
it has again dropped to less than 2 for the heat equation but about 190 for the fractional case. Of course, for 
T > T* the situation completely reverses. With J — 10, a = 1/4 and T = 10 T* the growth factor is a possibly 
workable value of around 600; while for the heat equation it is greater than 10 25 . We reiterate that the apparent 
contradiction between the theoretical ill-conditioning and numerical stability is due to the spectral cutoff present 
in any practical reconstruction procedure. 

Next we examine the influence of the fractional order a on the inversion step more closely. To this end, we 
expand the initial condition v in the piecewise linear finite element basis functions defined on a uniform partition 
of the domain Q = (0,1) with 100 grid points. Then we compute the discrete forward map F from the initial 
condition to the final time data g = u(T), defined on the same mesh. Numerically, this can be achieved by a fully 
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Table 1. The critical values T* for fractional backward diffusion. 


a\J 

3 

5 

10 

1/4 

0.0442 

0.0197 

0.0059 

1/2 

0.0387 

0.0163 

0.0049 

3/4 

0.0351 

0.0142 

0.0040 


discrete scheme based on the LI approximation in time and the finite difference method in space; see Appendix 
roi for a description of the numerical method. The ill-posed behavior of the discrete inverse problem is then 
analyzed using singular value decomposition. A similar experimental setup will be adopted for other examples 
below. 




k 

(b) singular value spectrum 


Figure 3. (a) The condition number v.s. the fractional order a, and (b) the singular value 
spectrum at T = 0.01 for the backward fractional diffusion. We only display the first 50 singular 
values. 


The numerical results are shown in Fig. [3] The condition number of the (discrete) forward map F stays mostly 
around O(10 4 ) for a fairly broad range of a values, which holds for all three different terminal times T. This can 
be attributed to the fact that for any a £ (0,1), backward fractional diffusion amounts to a two spacial derivative 
loss, cf. (13721) . Unsurprisingly, as the fractional order a approaches unity, the condition number eventually blows 
up, recovering the severely ill-posed nature of the classical backward diffusion problem, cf. Fig. [3]) a). Further, 
we observe that the smaller is the terminal time T, the quicker is the blowup. The precise mechanism for this 
observation remains unclear. Interestingly, the condition number is not monotone with respect to the fractional 
order a, for a fixed T. This might imply potential nonuniqueness in the simultaneous recovery of the fractional 
order a and the initial data v. The singular value spectra at T = 0.01 are shown in Fig. 0b). Even though 
the condition numbers for a — 1/4 and a = 1/2 are quite close, their singular value spectra actually differ 
by a multiplicative constant, but their decay rates are almost identical, thereby showing comparable condition 
numbers. This shift in singular value spectra can be explained by the local decay behavior of the Mittag-Leffler 
function, cf. Fig. 0a): the smaller is the fractional order a, the faster is the decay around t = 0. 

Even though the condition number is very informative about the (discretized) problem, it does not provide a full 
picture, especially when the condition number is large, and the singular value spectrum can be far more revealing. 
The spectra for two different a values are given in Fig. 0 At a = 1/2, the singular values decay at almost the 
same algebraic rate, irrespective of the terminal time T. This is expected from the two-derivative loss for any 
a < 1. However, for a = 1, the singular values decay exponentially, and the decay rate increases dramatically 
with the increase of the terminal time T. For T = 0.001, there are a handful of “significant” singular values, 
say above 10~ 3 , but when the time T increases to T = 1, there is only one meaningful singular value remaining. 
The distribution of the singular values has important practical consequences. For a small time T, the first few 
singular values for the classical diffusion case actually might be much larger than that for the fractional case, 
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which indicates that the classical case is actually numerically much easier to recover in this regime, concurring 
with the observations drawn from Table [l] For example, at T = 0.001, the first twenty singular values are larger 
than the fractional counterpart, cf. Fig. 0b), and hence, the first twenty modes, i.e., left singular vectors, are 
more stable in the reconstruction procedure. 




(a) a = 1/2 


(b) a = 1 


Figure 4. The singular value spectrum of the forward map F from the initial data to the final 
time data, for (a) a = 1/2 and (b) a = 1, at four different times for the time fractional backward 
diffusion. 


The mathematical model iim is rescaled with a unit diffusion coefficient. In practice, there is always a diffusion 
coefficient a in the elliptic operator, i.e., 

dtU — V • (crVit) + qu = f. 

For example, the thermal conductivity a of the gun steel at moderate temperature is about 1.8 x ICY 5 m 2 /s m 
and the diffusion coefficient of the oxygen in water at 25 Celsius is 2.10 x 10 J cm 2 /s [T8] . Mathematically, this 
does not change the ill-posed nature of the inverse problem. However, the presence of a diffusion coefficient a has 
important consequence: it enables the practical feasibility of the classical backward diffusion problem (and likely 
for many other inverse problems for the diffusion equation). Physically, a constant conductivity a amounts to 
rescaling the final time T by T' = T/cr. In the fractional case, a similar but nonlinear scaling law T' = T/a 1 ^ 01 
remains valid. 

Numerically, time fractional backward diffusion has been extensively studied. Liu and Yamamoto m proposed 
a numerical scheme for the one-dimensional fractional backward problem based on the quasi-reversibility method 
153 . and derived error estimates for the approximation, under a priori smoothness assumption on the initial 
condition. This represents one of the first works on inverse problems in anomalous diffusion. Later, Wang and 
Liu [99] studied total variation regularization for two-dimensional fractional backward diffusion, and analyzed its 
well-posedness of the optimization problem and the convergence of an iterative scheme of Bregman type. Wei 
and Wang m developed a modified quasi-boundary value method for the problem in a general domain, and 
established error estimates for both a priori and a posteriori parameter choice rules. In view of better stability 
results in the fractional case, one naturally expects better error estimates than the classical diffusion equation, 
which is confirmed by these studies. 

3.2. Sideways fractional diffusion. Next we consider the sideways problem for time fractional diffusion. There 
are several possible formulations, e.g., the quarter plane and the finite space domain. The quarter plane sideways 
fractional diffusion problem is as follows. Let u(x, t) be defined in (0, oo) x (0, oo) by 

d^u — u xx = 0, x > 0, t > 0, 

and the boundary and initial conditions 


u(a;,0)=0 and m(0, t) = f(t), 
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where we assume \u{x,t)\ < cie C2X . We do not know the left boundary condition /, but are able to measure u 
at an intermediate point x = L > 0, h(t) = u(L,t). The inverse problem is: given the (noisy) data h, find the 
boundary condition /. The solution u of the forward problem is given by a convolution integral with the kernel 
being the spatial derivative K a , x (x, s) of the fundamental solution K a (x, s), cf. (12.211 . by 

u{x,t) = / K a , x {x,t — s)f(s)ds. 

Jo 

This representation is well known for the case a = 1, and it was first derived by Carasso HU; see also [8] for 
related discussions. It leads to a convolution integral equation for the unknown / in terms of the given data h 

h(t) = [ R a (t — s)f(s)ds, 

Jo 

where the convolution kernel R a (s) is given by a Wright function in the form 

R ^ = ^ W ~ 

= y' (~ L ) fc s - fe f 
£jfc!T(-ffc + 2—f) 

In case of a = 1, i.e., classical diffusion, the kernel R(s) is given explicitly by 

R(s) = ^ s ~ z e~~ £ C°°( 0, oo). 

2 y 7T 

Since all its derivatives vanish at s = 0, the classical theory of Volterra integral equations of the first kind 
[56] implies the extreme ill-conditioning of the problem. This is not surprising: We are, after all, mapping a 
function / £ C°(0,oo) to an element in C°°(0, oo). The conditioning of the time fractional sideways problem 
again depends on the convolution kernel R a and its derivatives at s = 0 and in this case is the value of the Wright 
function W- a / 2 , 2 -o,/ 2 (—z) and its derivatives as z -> oo. These are again zero (in fact the Wright function 
W- a / 2 , 2 -a/ 2 (z ) also decays exponentially to zero for large negative arguments, cf. Lemma 12.211 . and thus the 
fractional sideways problem is also severely ill-posed. However, this analysis does not show their difference in the 
degree of ill-posedness: even though both are severely ill-posed, their practical computational behavior can still 
be quite different, as we shall see below. 

To see their difference in the degree of ill-posedness, we examine another variant of the sideways problem on 
a finite interval = (0,1), with Cauchy data prescribed on the axis x = 0, i.e. given zero initial condition 
uo = u(x, 0) = 0, recovering h = u(l, t) from the lateral Cauchy data at x = 0: 

u(0,t) = f(t), u x (0,t) = g(t), t> 0 


This problem is also known as the lateral Cauchy problem in the literature. In the case a = 1, it is known that 
the inverse problem is severely ill-posed mm- To gain insight into the fractional case, we apply the Laplace 
transform. With ~ being the Laplace transform in time, and noting the Laplace transform of the Caputo derivative 
C{dtU ) = z a u{z) — z a ~ 1 uo H3I Lemma 2.24], we deduce 


z a u(x,z)-u xx (x,z) = 0, n(0) =/, u x (0)=g. 

The general solution u is given by u(x,z) = f coshz a ^ 2 x + fl 51 ”!),/— and thus the solution h(z) = u(l,z) at 


x = 1 is given by 


h = f cosh z a/2 + 


sinh z a ' 2 
^ z“/ 2 


The solution h{t) can then be recovered by an inverse Laplace transform 

1 


(3.3) 


w)= ^L 


he zt dz, 


where Br = {z £ C : $tz = a, a > 0} is the Bromwich path. Upon deforming the contour suitably, this formula 
will allow developing an efficient numerical scheme for the sideways problem via quadrature rules EU, provided 
that the lateral Cauchy data is available for all t > 0. The expression (13731) indicates that, in the fractional case, 
the sideways problem still suffers from severe ill-posedness in theory, since the high frequency modes of the data 

a /2 

perturbation are amplified by an exponentially growing multiplier e 2 . However, numerically, the degree of 
ill-posedness decreases dramatically as the fractional order a decreases from unity to zero, since as a —> 0 + , the 
multipliers are growing at a much slower rate, and thus we have a better chance of recovering many more modes 
of the boundary data. In other words, both the classical and fractional sideways problems are severely ill-posed 
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in the sense of error estimates between the norms in the data and unknowns; but with a fixed frequency range, 
the behavior of the time fractional sideways problem can be much less ill-posed. Hence, anomalous diffusion 
mechanism does help substantially since much more effective reconstructions are possible in the fractional case. 



Figure 5. (a) The condition number and (b) singular value spectrum at T = 1 for the time 
fractional sideways problem. 


Next we illustrate the point numerically. The numerical results for the sideways problem are given in Fig. [5j 
It is observed that the degree of ill-posedness of the finite-dimensional discretized version of the inverse problem 
indeed decreases dramatically with the decrease of the fractional order a, cf. Fig. EJa), which agrees well with 
the preceding discussions. Surprisingly, for T = 1 there is a sudden transition around a = 1/2, below which the 
sideways problem behaves as if nearly well-posed, but above which the conditioning deteriorates dramatically with 
the increase of the fractional order a and eventually it recovers the properties of the classical sideways problem. 
Similar transitions are observed for other terminal times. This might be related to the discrete setting, for which 
there is an inherent frequency cutoff. Further, as the fractional order a approaches zero, the problem reaches a 
quasi-steady state much quicker and thus the forward map F can have only fairly localized elements along the 
main diagonal. To give a more complete picture, we examine the singular value spectrum in Fig. [5jb). Unlike the 
backward diffusion problem discussed earlier, the singular values are actually decaying only algebraically, even for 
a = 1, and then there might be a few tiny singular values contributing to the large condition number. The larger 
is the fractional order a, the more tiny singular values are in the spectrum. Hence, in the discrete setting, even 
for a = 3/4, the problem is still nearly well-posed, despite the large apparent condition number, since a few tiny 
singular values with a distinct gap from the rest of the spectrum are harmless in most regularization techniques. 

Physically this can also be observed in Fig. [6] where the forward map F is from the Dirichlet boundary 
condition x = 1 to the flux boundary condition at x = 0, in a piecewise linear finite element basis. Pictorially, the 
forward map F is only located in the upper left corner and has a triangular structure, which reflects the casual or 
Volterra nature of the sideways problem for the fractional diffusion equation. We note that the casual structure 
should be utilized in developing reconstruction techniques, via, e.g., Lavrentiev regularization 1561 . For small a 
values, e.g., a = 1/4, the finite element basis at the right end point x = 1 is almost instantly transported to the 
left end point x = 0, whose magnitude is slightly decreased, but with little diffusive effect, resulting a diagonally 
dominant forward map. However, as the fractional order a increases towards unity, the diffusive effect eventually 
kicks in, and the information spreads over the whole interval. Further, for large a values, it takes much longer 
time to reach the other side and there is a lag of information arrival, which explains the presence of tiny singular 
values. The larger is the fractional order a, the smaller is the magnitude, i.e., the less is the amount of the 
information reached the other side. Hence, one feasible approach is to recover only the boundary condition over a 
smaller subinterval of the measurement time interval. This idea underlies one popular engineering approach, the 
sequential function specification method @J [64]. 

The sideways problem for the classical diffusion has been extensively studied, and many efficient numerical 
methods have been developed and analyzed [111 f8i l24l 1251 . In the fractional case, however, there are only a 
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(a) a = 1/4 
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(b) a = 1/2 



0.4 

0.2 

0 

1 
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0.06 

0.04 

0.02 

0 

1 


(c) a = 3/4 (d) a = 1 

Figure 6. The Jacobian map F for a = 1/4, 1/2, 3/4 and 1, from the interval (0,1) itself. 


few works on numerical schemes, mostly for one-dimensional problems, and there seems no theoretical study 
on stability etc. Murio ESI EH developed several numerical schemes, e.g., based on space marching and finite 
difference method, for the sideways problem, but without any analysis. Qian [85] discussed about the ill-posedness 
of the quarter plane formulation of the sideways problem using the Fourier analysis, based on which a mollifier 
method was proposed, with error estimates provided. In [57], the recovery of a nonlinear boundary condition from 
the lateral Cauchy data was studied using an integral equation approach, and a convergent fixed point iteration 
method was suggested. The influence of the imprecise specification of the fractional order a on the reconstruction 
was examined. Zheng and Wei m proposed a mollification method for the quarter plane formulation of the 
sideways problem, by convoluting the fractional derivative with a smooth kernel, and derived error estimates for 
the approximation, under a prior bounds on the solution. The Cauchy problem of the time fractional diffusion 
has been numerically studied in Ena ■ In particular, with the separation of variables, a Volterra integral equation 
reformulation of the problem was derived, from which the ill-posedness of the Cauchy problem follows directly. 
All these works are concerned with the one-dimensional case, and the high dimensional case has not been studied. 

3.3. Inverse source problem. A third classical linear inverse problem for the diffusion equation is the inverse 
source problem, i.e., the recovery of the source term / from lateral boundary data or final time data. Clearly, one 
piece of boundary data or final time data alone is insufficient to uniquely determine a general source term, due to 
dimensionality disparity. To restore the possible uniqueness, as usual, we look for only a space- or time-dependent 
component of the source term /. With different combinations of the data and source term, we get several different 
(and not equivalent) formulations of the inverse source problems. Below we examine several of them briefly. By 
the linearity of the forward problem, we without loss of generality, assume a zero initial data v = 0 and a zero 
potential q — 0 throughout this part. 

First, suppose we can measure the solution u at the final time t = T, and aim at recovering either a space 
dependent or time dependent component of the source term /. Like before, we resort to the separation of variables. 
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For the case of a space dependent only source term /(x), the solution u to the forward problem is given by 


00 rt 

3=1 Jo 


){fAj)<t>3 dT 


= E v x (i - 

3 = 1 

Hence the measured data g = u(T) is given by 


E a , i(-A J T“))(/,^)^. 

A 7 


J=1 


By taking inner product with (f>j on both sides, we arrive at the following representation of the source term / in 
terms of the measured data g 


(3.4) 


/ = !> 


(s>&) 


f=i 


l-£ a ,i(-AjT Q 0 


By the complete monotonicity of the Mittag-Leffler function E a ^(—t) on the positive real axis R + [82], we deduce 

1 > £a,,i(-AiT a ) > £ Q ,i(-A 2 T a ), 

and thus the formula is well defined for any T > 0, and gives the precise condition for the existence of a 

source term. Even with a modest value of the terminal time T, the factor 1 — E a ,i(—AiT“) is close to unity for all 
small a values, especially for those close to zero. Each frequency component (/, (f>j) differs from (g, (f>j) essentially 
by a factor Aj, which amounts to two derivative loss in space. Actually one can show 


ll/lli 2 (n) < c ll5 , llir 2 (n)- 

This behavior is identical with that for the backward fractional diffusion. The statement holds also for the inverse 
source problem for the classical diffusion case. This is not surprising, since with a space dependent source term 
/, the solution u to the forward problem can be split into the steady solution u 3 and the decaying solution Ud, 
i.e., u = u s + Ud, where u 3 and Ud solve 

-u” = f, u s (0) = u s (l) = 0. 


and 

d?u d -u d ,xx = 0, Ud(0, x) = f(x), Ud(0,t) = Ud(l,t) = 0, 

respectively. By the decay behavior of the solution Ud, the steady state component u 3 is dominating, which 
amounts to a two spatial derivative loss. This is fully confirmed by the numerical experiments, cf. Fig. [3 It is 
observed that the condition number is almost independent of the fractional order a, and it is of order O(10 3 ), 
reflecting the mildly ill-posed nature of the inverse problem. In particular, for large terminal time T, the singular 
value spectra are almost identical for all fractional orders, decaying to zero at an algebraic rate, cf. Fig. 0b). 

Next we turn to the time dependent case, i.e., seeking a source term / of the form ,f(x,t) = p(t)q(x), with 
a known spacial component q(x), from the final time data g = u(T). Mathematically, the inverse problem 
even for the classical diffusion equation has not been completely analyzed. The inclusion of a nontrivial term 
q(x) is important since without this there is nonuniqueness. To see this, we take u to satisfy ut — u xx = f(t) 
on (0,1) x (0 ,T) with initial data w(®,0) = 1 and a homogeneous Neumann boundary condition —u x (0,t) = 
u x (l,t) = 0. Then one solution satisfying u(x,T) = g(x) = 1 is given by u(x,t) = 1 and / = 0, but another is 
u(x,t) = cos( 27 t t/T) and / = (— 2n/T) sin(2nt/T). Likewise, in the fractional case, we can take u = cos(2ivt/T) 
for the second solution and define / to be its ath order Djrbashian-Caputo fractional derivative in time. 

Like previously, the solution u to m is given by 

CO pf 

U ( t ) = Yl / (* - ' r )“ _1 -Ec«,a(-Aj(t - T) a )p{T) d T{q,(t>j)<l>j. 

3 = 1 J ° 

Hence the measured data g(x) = u[x, T ) is given by 
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k 

(b) singular value spectrum 


Figure 7. Numerical results for the inverse source problem with final time data and a space 
dependent source term, (a) The condition number v.s. the fractional order a, and (b) singular 
value spectrum at T = 1. 




(a) T = 0.01 (b) T = 1 


Figure 8. The singular value spectrum at two different terminal times for the inverse source 
problem with a final time data at terminal time T, and f(x,t) = xp(t), an unknown time 
dependent component pit). 


By taking inner product with <p : j on both sides, we deduce 

(9,<f>i) = {<lAo) [ ( T-T) a ~ 1 E aia (-\ j {T-r) a )p(T)dr. 

Jo 

In the case of a = 1, the formula recovers the relation 

{gAj) = f e~ Xj(T - T) p{r)dT, 

Jo 

which resembles a finite-time Laplace transform or moment problem, and thus severely smoothing, which renders 
the inverse source problem severely ill-posed. Intuitively, the term can only pick up the information 

for t close to the terminal time T, and for t away from T, the information is severely damped, especially for high 
frequency modes, which leads to the severely ill-posed nature of the inverse problem. In the fractional case, the 
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forward map F from the unknown to the data is clearly compact, and thus the problem is still ill-posed. However, 
the kernel t 01 ^ 1 E a , a (—Xjt a ) is less smooth and decays much slower, and one might expect that the problem is less 
ill-posed than the canonical diffusion counterpart. To examine the point, we present the numerical results for the 
inverse problem in Fig. [8j It is severely ill-posed irrespective of the fractional order a: the singular values decay 
exponentially to zero without a distinct gap in the spectrum. In particular, for the terminal time T = 1, the 
spectrum is almost identical for all fractional orders a. For small T, the singular values still decay exponentially, 
but the rate is different: the smaller is the fractional order a, the faster is the decay, cf. Fig. Eta). Consequently, 
a few more modes of the source term p(r) might be recovered. In other words, due to a slower local decay of 
the exponential function e _At , compared with the Mittag-Leffler function t a ~ 1 E ai>a {— Af“), cf. Fig. QJa), actually 
more frequency modes can be picked up by normal diffusion than the fractional counterpart, cf. Fig. [8ta). This 
indicates that with sufficiently accurate data, at a small time instance, the sideways problem for normal diffusion 
may allow recovering more modes, i.e., anomalous diffusion does not help solve the inverse problem. 

In practice, the accessible data can also be the flux data at the end point, e.g., x = 0 or x = 1. We briefly 
discuss the case of recovering a time dependent component p(t) in the source term / = q(x)p(t) from the flux 
data at x = 0. By repeating the preceding argument, the data g := —u x (0,t) is related to the unknown p(t) by 

oo 

j=i Jo 

In 1881 Theorem 4.4], a stability result was established for the recovery of the time dependent component p(t). 
Along the same line of thought, under reasonable assumptions, one can deduce that 

l|p||c[0,T] < c\\dt g\\c[0,T]- 

The inverse problem roughly amounts to taking the ath order Djrbashian-Caputo fractional derivative in time. 
Hence as the fractional order a decreases from unity to zero, it becomes less and less ill-posed. For a close to 
zero, it is nearly well-posed, at least numerically. In other words, anomalous diffusion can mitigate the degree of 
ill-posedness for the inverse problem. To illustrate the discussion, we present in Fig. [9] some numerical results, 
where the forward map F is from the time dependent component p(t) to the flux data g{t) at x — 0, both defined 
over the interval [0, T], discretized using a continuous piecewise linear finite element basis. The condition number 
of the discrete forward map F decreases monotonically as the fractional order a decreases from unity to zero, 
confirming the preceding discussions. Further, the terminal time T does not affect the condition number to a 
large extent. 




(b) singular values 


Figure 9. Numerical results for the inverse source problem with flux data at x = 0 and f(x, t) = 
xp(t), an unknown time dependent component p(t). (a) The condition number of the discrete 
forward map and (b) singular value spectrum at T = 1. 


It is widely accepted in inverse heat conduction that an inverse problem will be severely ill-posed when the 
data and unknown are not aligned in the same space/time direction, and only mildly ill-posed when they do align 
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with each other. Our discussions with the inverse source problems indicate that the observation remains valid in 
the time fractional diffusion case. In particular, although not presented, we note that the inverse source problem 
of recovering a space dependent component from the lateral Cauchy data is severely ill-posed for both fractional 
and normal diffusion. In the simplest case of a space dependent only source term, it is mathematically equivalent 
to unique continuation, a well known example of severely ill-posed inverse problems. 

The inverse source problems for the classical diffusion equation have been extensively studied; see e.g., [3 E3 ID- 
Inverse source problems for FDEs have also been numerically studied. Zhang and Xu EH established the unique 
recovery of a space dependent source term in EH) with pure Neumann boundary data and overspecified Dirichlet 
data at x = 0. This is achieved by an eigenexpansion and Laplace transform, and the uniqueness follows from a 
unique continuation principle of analytic functions. Sakamoto and Yamamoto [89] discussed the inverse problem 
of determining a spatially varying function of the source term by final overdetermined data in multi-dimension, 
and established its well-posedness in the Hadamard sense except for a discrete set of values of the diffusion 
constant, using an analytic Fredholm theory. Very recently, Luchko et al [68] showed the uniqueness of recovering 
a nonlinear source term from the boundary measurement, and developed a numerical scheme of fixed point 
iteration type. Aleroev et al [2] showed the uniqueness of recovering a space dependent source term from integral 
type observational data. Recently, there are many numerical studies on this class of inverse problems. In m, 
the numerical recovery of a spatially varying function of the source term from the final time data in a general 
domain was studied using a quasi-boundary value problem method; see also mm related studies. Wang et al 
[TOO | proposed to determine the space-dependent source term from the final time data in multi-dimension using 
a reproducing kernel Hilbert space method. 

3.4. Inverse potential problem. Now we consider a nonlinear inverse coefficient problem for the time fractional 
diffusion equation: given the final time data g = u(T), find the potential q in the model 

(3.5) d?u — u xx + qu = 0 in Q, 

with a homogeneous Neumann boundary condition and initial data v. The parabolic counterpart has been 
extensively studied 13811151 [16] . where it was shown that the problem is nearly well-posed in the Hardamard sense 
in suitable Holder space, under certain conditions, using the strong maximum principle. In [35], an elegant fixed 
point method was developed, and the monotone convergence of the method was established. It can be adapted 
straightforwardly to the fractional case: given an initial guess q° , compute the update q k recursively by 

jb+i _ 9" - dtu(x, T ; q k ) 

Q 

9 

where the notation u(x,T-q k ) denote the solution to problem (1331) with the potential q k at t = T. Since the 
strong maximum principle is still valid for the time fractional diffusion equation m , the scheme is monotonically 
convergent, under suitable conditions. 

As the terminal time T —» oo, the problem recovers a steady-state problem, and the scheme amounts to twice 
numerical differentiation in space and converges within one iteration, provided that the data g is accurate enough. 
Hence, it is natural to expect that the convergence of the scheme will depends crucially on the time T : the larger 
is the time T, the closer is the solution u to the steady state solution; and thus the faster is the convergence of 
the fixed point scheme. By Lemma 12.11 as the fractional order a approaches zero, the solution u decays much 
faster around t = 0 than the classical one, i.e., a = 1. In other words, the fractional diffusion problem can reach 
a “quasi-steady state” much faster than the classical one, especially for a close to zero, and the scheme will then 
converge much faster. 

To illustrate the point, we present in Fig. [I0]some numerical results of reconstructing a discontinuous potential 
q = l + 2a;X[o,o.5] +2(1 — Xto.s.i]*) (with \s being the characteristic function of the set S) from exact data, in order 
to illustrate the convergence behavior of the fixed point scheme. In the figure, e denotes the relative L 2 (Q) error. 
The numerical results fully confirm the preceding discussions: at a fixed time T, the smaller is the fractional 
order a, the faster is the convergence; and at fixed a , the larger is the time T, the faster is the convergence. 
Numerically, one also observes the monotone convergence of the scheme. 

Generally, the recovery of a coefficient in FDEs has not been extensively studied. Cheng et al m established 
the unique recovery of the fractional order a G (0,1) and the diffusion coefficient from the lateral boundary 
measurements. It represents one of the first mathematical works on invere problems for FDEs, and has inspired 
many follow-up works. Yamamoto and Zhang EH established conditional stability in determining a zeroth-order 
coefficient in a one-dimensional FDE with one half order Caputo derivative by a Carleman estimate. Carleman 
estimates for time fractional diffusion were discussed in (10811131162] . In |73| . the unique determination of the 
spatial coefficient and/or the fractional order from the data on a subdomain was shown for a positive initial 
condition. Wang and Wu [97] studied the simultaneous recovery of two time varying coefficients, i.e., a kernel 
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(a) error at T = 0.1 (b) error at a = 1/2 

Figure 10. Numerical results, i.e., the relative L 2 (Q) error e, for the inverse potential problem 
from exact final time data at (a) T = 0.1 and (b) a = 1/2. 


function and a source function, from the additional integral observation in multi-dimension, using a fixed point 
theorem. All these works are concerned with the theoretical analysis, ant there are even fewer works on the 
numerical analysis of related inverse problems. Li et al j55] suggested an optimal perturbation algorithm for 
the simultaneous numerical recovery of the diffusion coefficient and fractional order in a one-dimensional time 
fractional FDE. In [50], the authors considered the identification of a potentiai term from the iateral flux data at 
one fixed time instance corresponding to a complete set of source terms, and established the unique determination 
for “small” potentials. Further, a Newton type method was proposed in [50], and its convergence was shown. 

Even though our discussions have focused on time fractional diffusion, which invoives one single fractional 
derivative in time, it is also possible to consider equations where the time derivative involves multiple factional 
orders, i.e., Ylk=i for a sequence cti > «2 > ... > a m [Ml 2D; see ES f° r some first uniqueness results 

for inverse coefficient probiems in the multi-dimensional case. Further extensions include the distributed-order, 
spatially and/or temporally variable-order and tempered fractional diffusion, to better capture certain physicai 
processes, for which, however, reiated inverse problems have not been discussed at all. 

3.5. Fractional derivative as an inverse solution. One of the very first undetermined coefficient probiems 
for PDEs was discussed in the paper by B. Frank Jones Jr. [52] (see also [8] Chapter 13]). This is to determine 
the coefficient a(t) from 

ut = a(t)u xx , 0 < x < oo, t > 0 
u(a 5,0) = 0, — a(t)u x (0, t) = g(t), 0 <t<T 

under the over-posed condition of measuring the “temperature” at x = 0 

u(0,£) = ip(t) 

In [52] . B. Frank Jones Jr. provided a complete analysis of the problem, by giving necessary and sufficient 
conditions for a unique soiution as well as determining the exact level of ill-conditioning. The key step in the 
analysis is a change of variables and conversion of the problem to an equivalent integral equation formulation. 
Perhaps surprisingly, this approach involves the use of a fractional derivative as we now show. 

The assumptions are that g is continuous and positive and ip is continuously differentiable with ip( 0) = 0 and 
i)' > 0 on (0, T). In addition, the function h(t) defined by 

h( t) = __ 

fo (* - 

satisfies limt_>o h(t) = h o > 0. Note that h is the ratio of the two data functions; the flux g and the Djrbashian- 
Caputo derivative of order 1/2 of ip. If we define hi = inf h and h s = sup h on [0, T] and look at the space 
G := {a G C[0, T) : h 2 < a(t) < h 2 }, then it was shown that any a G G satisfies the inverse problem must also 
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solve the integral equation 


a(t) 


__ 

fo a(s) ds]- 1 / 2 dr 


= : Ta, 


and vice-versa. The main result in [52] is that the operator T has a unique fixed point on Q and indeed T is 
monotone in the sense of preserving the partial order on Q, i.e., if a\ > 02 then Tai < Taz- 

Given these developments, it might seem that a parallel construction for the time fractional diffusion coun¬ 
terpart, dtU — a(t)u X x, would be relatively straightforward but this seems not to be the case. The basic steps 
for the parabolic version require items that just are not true in the fractional case, such as the product rule, and 
without these the above structure cannot be replicated or at least not without some further ingenuity. 


4. Inverse problems for space fractional diffusion 

Now we turn to differential equations involving a fractional derivative in space. There are several possible 
choices of a fractional derivative in space, e.g., Djrbashian-Caputo fractional derivative, Riemann-Liouville frac¬ 
tional derivative, Riesz derivative, and fractional Laplacian [3], which all have received considerable attention. In 
recent years, the use of the fractional Laplacian is especially popular in high-dimensional spaces, and admits a 
well-developed analytical theory. We shall focus on the left-sided Djrbashian-Caputo fractional derivative g'D.f, 
P £ (1, 2), and the one-dimensional case, and consider the following four inverse problems: inverse Sturm-Liouville 
problem, Cauchy problem for a fractional elliptic equation, backwards diffusion, and sideways problem. 

4.1. Inverse Sturm-Liouville problem. First we consider the following Sturm-Liouville problem on the unit 
interval SI = (0, 1): find u £ Hq(Q) n H^( SI) and A £ C such that 

(4.1) — oDxU + qu = \u in SI, 

with a homogeneous Dirichlet boundary condition u(0) = u(l) = 0. A Sturm-Liouville problem of this form 
was considered by Mkhitar M Djrbashian [22l ll9l in 1960s to construct certain biorthogonal basis for spaces of 
analytic functions; see also [78]. Like before, with /3 = 2, it recovers the classical Sturm-Liouville problem. In 
the case of a general potential q, in the fractional case, little is known about the analytical properties of the 
eigenvalues and eigenfunctions. For the case of a zero potential q — 0, there are countably many eigenvalues {Aj} 
to (14.11) . which are zeros of the Mittag-Leffler function E f s, 2 (—A). The corresponding eigenfunctions are given by 
xEp^i—XjX 13 ). Using the exponential asymptotics on the Mittag-Leffler function in Lemma [2.11 one 1931150) can 
show that asymptotically, the eigenvalues A, are distributed as 

I'M ~ (27 vj) p and arg(Aj) - ^ ^ . 

Hence, for any /3 £ (1,2), there are only a finite number of real eigenvalues to (14.11) . and the rest appears as 
complex conjugate pairs. 

It is well known that eigenvalues contain valuable information about the boundary value problem. For example 
it is known that the sequence of Dirichlet eigenvalues can uniquely determine a potential q symmetric with respect 
to the point x = 1/2, and together with additional spectral information, one can uniquely determine a general 
potential q\ see I121I86I for an overview of results on the classical inverse Sturm-Liouville problem. In the fractional 
case, the eigenvalues are generally genuinely complex, and a complex number may carry more information than 
a real one. Thus one naturally wonders whether these complex eigenvalues do contain more information about 
the potential. Numerically the answer is affirmative. To illustrate this, we show some numerical reconstructions 
in Fig. 1111 obtained by using a frozen Newton method and representing the sought-for potential q in Fourier 
series [50;. The Dirichlet eigenvalues can be computed efficiently using a Galerkin finite element method l45l . 
One observes that one single Dirichlet spectrum can uniquely determine a general potential q. Unsurprisingly, as 
the fractional order /? tends two, the reconstruction becomes less and less accurate, since in the limit /3 = 2, the 
Dirichlet spectrum cannot uniquely determine a general potential q. Theoretically, the surprising uniqueness in 
the fractional case remains to be established. 

Naturally, one can also consider the Riemann-Liouville case: 

(4.2) —^D^u + qu = Au in SI, 

with u(0) = u(l) = 0. Like before, little is known about the analytical properties of the eigenvalues and 
eigenfunctions. In the case of a zero potential q = 0, there are countably many eigenvalues to (1X21) . which are zeros 
of the Mittag-Leffler function Ep,p(— A), and the corresponding eigenfunctions are given by x^~ 1 Ep l p(—\jX^). 
Further, the asymptotics of the eigenvalues are still valid. Hence, for any ft £ (1, 2), there are only a finite number 
of real eigenvalues to (1X21) . and the rest appears as complex conjugate pairs. 
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(a) p = 5/3 (b) p = 7/4 

Figure 11. Numerical results for the inverse Sturm-Liouville problem with a Djrbashian- 
Caputo derivative for (a) P = 5/3 and (b) P = 7/4. The reconstructions are computed from the 
first eight eigenvalues (in absolute value) using a frozen Newton method [50i. 


The numerical results from the Dirichlet spectrum in the Riemann-Liouville case are shown in Fig. 1121 For a 
general potential q, the reconstruction represents only the symmetric part, which is drastically different from the 
Djrbashian-Caputo case, but identical with that for the classical Sturm-Liouville problem. Further, if we assume 
that the potential q is known on the left half interval, then the Dirichlet spectrum allows uniquely reconstructing 
the potential q on the remaining half interval, cf., Fig. 1121 b). These results indicate that in the Riemann-Liouville 
case the complex spectrum is not more informative than the classical Sturm-Liouville problem. The precise 
mechanism underlying the fundamental differences between the Djrbashian-Caputo and Riemann-Liouville cases 
awaits further study. However, as in the classical case P = 2, one can show that the linearized derivative of the 
map q —> u( 1; A, q) around q = 0 cannot span more than the subspace of even functions in L 2 ( fl). 




(a) whole interval (b) half interval 


Figure 12. Numerical results for the inverse Sturm-Liouville problem with a Riemann-Liouville 
fractional derivative of order P = 4/3. The reconstructions are computed from the first eight 
eigenvalues (in absolute value) using a frozen Newton method [50] . 
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In general, the Sturm-Liouville problem with a fractional derivative remains completely elusive, and numerical 
methods such as finite element method [46] provide a valuable (and often the only) tool for studying its analyt¬ 
ical properties. For a variant of the fractional Sturm-Liouville problem, which contains a fractional derivative 
in the lower-order term, Malamud m established the existence of a similarity transformation, analogous to the 
well-known Gel’fand-Levitan-Marchenko transformation, and also the unique recovery of the potential from mul¬ 
tiple spectra. In the classical case, the GeLfand-Levitan-Marchenko transformation lends itself to a constructive 
algorithm [86]; however, it is unclear whether this is true in the fractional case. In [50], the authors proposed a 
Newton type method for reconstructing the potential, which numerically exhibits very good convergence behavior. 
However, a rigorous convergence analysis of the scheme is still missing. Further, the uniqueness and nonuniqueness 
issues of related inverse Sturm-Liouville problems are outstanding. Last, as noted above, there are other possible 
choices of the space fractional derivative, e.g., fractional Laplacian and Riesz derivative. It is unknown whether 
the preceding observations are valid for these alternative derivatives. 

4.2. Cauchy problem for fractional elliptic equation. One classical elliptic inverse problem is the Cauchy 
problem for the Laplace equation, which plays a fundamental role in the study of many elliptic inverse problems 
w ■ A first example was given by Jacques Hadamard m to illustrate the severe ill-posedness of the Cauchy 
problem, which motivated him to introduce the concept of well-posedness and ill-posedness for problems in 
mathematical physics. So a natural question is whether the Cauchy problem for the fractional elliptic equation 
is also as ill-posed? To illustrate this, we consider the following fractional elliptic problem on the rectangular 
domain fi = {(a:, y) £ R 2 : 0 < x < 1, 0 < y < 1} 

(4.3) aQfu + oDyU = 0 in fl, 

with the fractional order 0 £ (1,2) and the Cauchy data 

du 

u(x, 0) = g(x) and — (x,0) = h(x), 0 < x < 1 

where v is the unit outward normal direction. With 0 = 2, it recovers the Cauchy problem for the Laplace 
equation. By applying the separation of variables, we assume that u(x,y) = cp(x)ip(y), which directly gives for 
some scalar A £ C that 

o D £<p(x) = —X(p(x), 
oDy ip(y) = A ip{y), 

Let (A j,4>j) be a Dirichlet eigenpair of the Caputo derivative operator — on the unit interval D = (0,1), i.e., 

< pj(x ) = xEp t 2 (—AjO^), and |Ay| —I oo as j —>• oo | 51 |: see Section [47T] for further details. With the choice cp = <f>j 
and the Cauchy data pair ( g , hj) = (0, —xEp^{XjX^)/Xj), the component ipj satisfies 

oDy'tpj(y) = A ji>j(y) in y £ (0,oo), 

with the initial condition ipj (0) = 0 and -^ipj{ 0) = 1/A j. Using the relation -j^x 1 ^ 1 Ep tl {Xx^) = Xx 1 ^ 2 Ep^-i(Xx^) 
[531 pp. 46], we deduce that the solution ipj to the fractional ordinary differential equation is given by 

fpjiv) = yEpAXjy p )/Xj. 

Hence, Uj[x,y ) = xEp t 2 (—XjX l3 )yEp^(Xjy l3 )/X 2 is a solution to the Cauchy problem with g = 0 and hj(x) = 
—xEp t 2 {—XjX^)/Xj. By the exponential asymptotics of the Mittag-Leffler function, cf. Lemma 12.11 we deduce 
that hj(x) —I 0 as j —>• oo, whereas for any y > 0, the solution Uj(x,y) oo as j oo, in view of the exponential 
growth of the Mittag-Leffler function Ep^{z), cf. Lemma 12.11 This indicates that the Cauchy problem for the 
fractional elliptic equation is also exponentially ill-posed. However, the interesting question of the degree of ill- 
posedness, in comparison with the classical case, is unclear and certainly worthy of further study. Further, we 
note that the numerical solution of the fractional elliptic equation m is highly nontrivial, and there seems no 
efficient yet rigorous solver available in the literature, due to a lack of solution theory for such problems. 

4.3. Backward problem. Now we return to the backward diffusion problem with fractional derivatives in the 
space variable(s). Let Q = (0,1) be the unit interval. Then the one-dimensional space fractional diffusion equation 
is given by 

u t — oD^u = 0, (x,t) £ B x (0, oo), 

where the fractional order 0 £ (1, 2). The equation is equipped with the following initial condition u(x, 0) = v and 
zero boundary condition u(0, t) = u(l, t) = 0. The backward problem is: given the final time data g(x) = u(x, T), 
find the initial data v. Since the Djrbashian-Caputo derivative operator qD£ with the zero Dirichlet boundary is 
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sectorial on suitable spaces |421 , the existence of a solution u follows from the analytic semigroup theory |8l[ 143] , 
and formally it can be represented by 

u(t) = e~ At v, 

where A is the representation of the Djrbashian-Caputo derivative operator — q on its domain. Formally, the 
solution v to the space fractional backward problem is given by 

AT 

v = e g. 

In case of /3 = 2, using the eigenpairs {(A-,, <frj)} and the L 2 (fl) orthogonality of the eigenfunctions it recovers 
the well known formula 

OO 

3 = 1 

The growth factor e XjT explains the severely ill-posed nature of the inverse problem. In the fractional case, such 
an explicit representation is no longer available since the corresponding eigenfunctions {4>j} are not orthogonal 
in L 2 (Q, ) (actually they can be almost linearly dependent), due to the non self adjoint nature of the Djrbashian- 
Caputo derivative operator — o'D.f ■ Nonetheless, according to the discussions in Section rm the eigenvalues {Aj} 
increase to infinity with the index j. and asymptotically lies on two rays. Hence, one naturally expects that the 
backward problem is also exponentially ill-posed. However, the magnitudes (and the real parts) of the eigenvalues 
grow at a rate slower than that of the standard Sturm-Liouvillc problem, and thus the space fractional backward 
problem is less ill-posed than the classical one. To illustrate the point, we present the numerical results in Fig. 
M For all fractional orders /3, the singular values decay exponentially, but the decay rate increases dramatically 
with the increase of the fractional order /? and the terminal time T. Hence, anomalous superdiffusion does not 
change the exponentially ill-posed nature of the backward problem, but numerically it does enable recovering 
more Fourier modes of the initial data v. 




(a) T = 0.01 (b) T = 0.1 

Figure 13. Numerical results for the space fractional backward diffusion problem, the singular 
value spectrum at two different time instances, (a) T = 0.01 and (b) T = 0.1. 


Last, we note that for other choices of the fractional derivative, e.g., the Riemann-Liouville fractional derivative 
and the fractional Laplacian |S|55], the magnitude of eigenvalues of the operator also tends to infinity, and the 
growth rate increases with the fractional order /3. Therefore, the preceding observations on the space fractional 
backward problem are expected to be valid for these choices as well. 

4.4. Sideways problem. Last we return to the classical sideways diffusion problem but now with a fractional 
derivative in space rather than in time. Let fl = (0,1) be the unit interval. Then the one-dimensional space 
fractional diffusion equation is given by 

Ut — qD^u = 0, ( x,t ) £ Q. x (0, oo), 
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where the fractional order (3 £ (1,2). The equation is equipped with an initial condition u(x, 0) = 0 and the 
following lateral Cauchy boundary conditions 

u(0, t) = f(t) and u x (0,t) = g(t), t > 0. 

We wish to compute the solution at x = 1, i.e., h(t) := u(l,t). In the case (3 — 2, the model recovers the 
standard diffusion equation, and we have already discussed the severe ill-conditioning of the classical case. Due to 
the nonlocal nature of the fractional derivative, one might expect that in the space fractional case, the sideways 
problem is less ill-posed. To see this, we take Laplace transform in time to arrive at (with ~ denoting the Laplace 
transform) 

zu(x, z ) - qD£u(x, z) = 0, 

with the initial conditions (at x = 0) 

u(0, z) = f(z) and u x ( 0, z) = g(z). 

The solution u(x, z) to the initial value problem is given by 

u(x,z) = f(z)Ep }1 (zx 0 ) +g{z)xEp, 2 (zx p ) 

and thus 

Hz) = f(z)Ef} : i(z) + g(z)Ef}' 2 (z). 

Like before, the boundary condition h(t) at x = 1 can be found by an inverse Laplace transform 

h(t) = f e zt h(z)dz. 

J Br 

In case of /3 = 2, this gives cosh y/z and sinh yfz/yfz multipliers to the data f(z) and g(z) resulting in the 
exponential ill-conditioning of the sideways heat problem. In the case of a general (3 £ (1,2), the exponential 
asymptotics in Lemma 12.II indicates that the problem still suffers from exponentially growing multipliers to the 
data, and thus the problem is still severely ill-conditioned. Simple computation shows that the multiplier is 
asymptotically larger for the fractional order /3 closer to unity. In other words, anomalous diffusion in space does 
not mitigate the ill-conditioned nature of the sideways problem, but actually worsens the conditioning severely. 

To further illustrate the point, we compute the forward map F from the Dirichlet boundary condition at x = 1 
to the flux at x = 0 numerically with a finite element in space/finite difference in time scheme, cf. Appendix IA.31 
for the details. The numerical results are presented in Fig. 1141 The singular value spectra clearly show the ill- 
posedness nature of the space fractional sideways problem: as the fractional order /3 increases from one to two, the 
majority of the singular values move upward, the decay of the singular values slows down, and thus the sideways 
problem becomes less and less ill-posed (but still severely so). Further, there are more tiny singular values kicking 
in as the fractional order (3 decreases to one, which indicates the inherent rank deficiency of the forward map 
F and might be relevant in the uniqueness of the inverse problem. This confirms the preceding analysis: the 
degree of ill-posedness worsens with the decrease of the fractional order /3, and the fractional counterpart is more 
ill-posed than the classical one. In other words, anomalous diffusion actually severely worsens the conditioning of 
the already very ill-posed sideways problem. Further, the numerical results tend to indicate that the Djrbashian- 
Caputo derivative with an order /3 £ (1,2) acts as an interpolation between the diffusion and convection, which 
results in a history mechanism in space: when the history piece runs from the left to the right, it is unlikely to 
propagate the information in the reverse direction; and the closer is the fractional order f3 to unity, the stronger 
is the directional effect. The latter is not counterintuitive, since in the limit of (3 = 1, the Djrbashian-Caputo 
fractional derivative qD^u recovers the first order derivative and the problem is of convection type, and surely 
no information can be convected backwards! 

In the case (3 = 2, one may equally measure the lateral Cauchy data at x = 1, and aims at recovering the 
Dirichlet boundary condition at x — 0. Clearly, this does not change the nature of the inverse problem, and it is 
equally ill-posed. Due to the directional nature of the Djrbashian-Caputo derivative o'D.f, one naturally wonders 
whether this “directional” feature does influence the ill-posed nature of the sideways problem. To illustrate the 
point, we repeat the preceding arguments and deduce 

zu(x, z) — qD£u(x, z) = 0, 

with the boundary conditions at x = 1 

u(l,z) = f{z) and u x (l, z) = g(z). 
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Figure 14. Singular value spectrum of the forward map F at times T = 0.1 and T — 1, for the 
sideways problem with Cauchy data at x = 0. 


To derive the solution, denote the initial conditions at x = 0 by f(z) = u(0,z) and g(z ) = u x (0,z). Then the 
solution u(x, z) to the initial value problem is given by 

u(x, z) = f{z)Epp{zx p ) + gxEp }2 (zx p ). 

Use the differentiation formula Ep : ~,(zx p ) = zx 1-2 Ep n _\(zx p ) [22 page 46], we deduce that at x = 1, 

there hold 

f{z)Ep i i(z) + gE Pt2 (z) = f(z), 
f{z)Epfi{z) + gEpp(z) = z~ 1 g(z). 

Solving the linear system yields the solution to the sideways problem 

j,' = fffl.i(-)/(") - z~ lE nAz)g{z) 

> Epp{zf - Ep }0 (z)Ep, 2 (z) ’ 


and accordingly the solution h(t) = u(0, t) is given by an inverse Laplace transform. The growth factors of the 
data / and g are Epp(z)/(Epp(z) 2 — Ep : o(z)Ep, 2 (z)) and z _1 Ep t2 (z)/(Epp(z) 2 — Ep : o(z)Ep, 2 (z)), respectively. 
The growth of these factors at large s argument determines the degree of ill-conditioning of the sideways problem. 
To this end, we appeal to the exponential asymptotic of the Mittag-Leffler function E a: p(z ), cf. Lemma [2.1l and 
note that Bromwhich path lies in the sector | argz| < n/2 to deduce that for large \z\, there holds 


Ep^zf 

Ep : o(z)Ep i2 (z) 


1 2 z 1 /? 2 Z l/H 

pr(i -p)z e ’ 

1 2 1 1 // 3-1 z 1 /? 

/3r(2-/3) 


Hence, the numerator Epp(z) 2 — Ep,o(z)Ep t2 (z) behaves like 

Epp(zf - Ep t0 (z)Ep t2 (z) - p T (2- p) zl,P ~ leZ ^ aS ^ ^ °°- 

This together with the exponential asymptotic of Epp(z) and Ep :2 (z) from Lemma |2. II indicates that the multi¬ 
pliers for / and ~g are growing at most at a very low-order polynomial rate, for large z. Hence, the high-frequency 
components of the data noise are not amplified much (at most polynomially instead of exponentially). The anal¬ 
ysis indicates that the sideways problem with the lateral Cauchy data specified at the point x = 1 is nearly 
well-posed, as long as the fractional order /3 is away from two, for which it recovers the classical ill-posed sideways 
problem for the heat equation. 

Next we illustrate the preceding discussions numerically. The behavior of the forward map F from the Dirichlet 
boundary at x = 0 to the flux data at x = 1 is shown in Fig. 1151 For a wide range of values of the fractional 
order /3, the condition number of the forward map F is of order 100, which is fairly mild, in view of the size of the 
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(a) condition number (b) T = 1 

Figure 15. Numerical results for the space fractional sideways problem, with the lateral Cauchy 
data at the point x = 1: (a) the condition number v.s. the fractional order /3 and (b) the singular 
value spectrum at T = 1. 


linear system, i.e., 100 x 100. When the fractional order /? increases towards two, the inverse problem recovers the 
classical sideways problem, and as expected, the condition number increases dramatically. However, the onset of 
the blowup depends on the terminal time T: the smaller is the time T, the smaller seems the onset value. The 
precise mechanism for this phenomenon is still unknown. For /3 < 7/4, the singular value spectrum only spans 
a narrow interval, resulting in a very small condition number. Physically, like before, this can be explained as 
the “convective” nature of the Djrbashian-Caputo fractional derivative: as the fractional order /3 tends to unity, 
the information at x = 0 is transported to x = 1, free from distortion, and thus the inverse problem is almost 
well-posed. In summary, depending on the location of the over-specified data, anomalous superdiffusion can either 
help or aggravate the conditioning of the sideways problem. 

Last, we would like to note that the study of space fractional inverse problems, either theoretical or numerical, 
is fairly scarce. This is partly attributed to the relatively poor understanding of forward problems for FDEs 
with a space fractional derivative: There are only a few mathematical studies on one-dimensional space fractional 
diffusion, and no mathematical study on multi-dimensional problems involving space fractional derivatives (of 
either Riemann-Liouville or Caputo type). Nonetheless, our preliminary numerical experiments show distinct 
new features for related inverse problems, which motivate their analytical studies. 

5. Concluding remarks 

Anomalous diffusion processes arise in many disciplines, and the physics behind is very different from normal 
diffusion. The unusual physics greatly influences the behavior of related forward problems. Further, it is well 
known that backward fractional diffusion is much less ill-posed than the classical backward diffusion, which has 
contributed to the belief that inverse problems for anomalous diffusion are always better behaved than that 
for the normal diffusion. In this work we have examined several exemplary inverse problems for anomalous 
diffusion processes in a numerical and semi-analytical manner. These include the sideways problem, backward 
problem, inverse source problem, inverse Sturm-Liouville problem and Cauchy problems. Our findings indicate 
that anomalous diffusion can give rise to very unusual new features, but they only partially confirm the belief: 
depending on the data and unknown, it may influence either positively or negatively the degree of ill-posedness 
of the inverse problem. 

The mathematical study of inverse problems in anomalous diffusion is still in its infancy. There are only a few 
rigorous theoretical results on the uniqueness, existence and stability, which mostly focus on the one-dimensional 
case, and there are many more open problems awaiting investigations. The development of stable and efficient 
reconstruction procedures is an active ongoing research topic. However, due to the nonlocality of the forward 
model, the construction of efficient schemes and their rigorous numerical analysis remain very challenging. This 
is especially true for space fractional FDEs, and there are almost no theoretical or rigorous numerical studies. 
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Appendix A. Numerical methods for special functions and FDEs 

A.l. Computation of the Mittag-Leffler and Wright functions. Like many special functions, the efficient 
and accurate numerical computation of the Mittag-Leffler function E at p(z) is delicate 1281 [2'9l [94l . An efficient 
algorithm relies on partitioning the complex plane C into different regions, where different approximations, i.e., 
power series, integral representation and exponential asymptotic for small values of the argument, intermediate 
values and large values, respectively, are used for efficient numerical computation; see [94] for the some partition 
and error estimates. The special case of the Mittag-Leffler function E ai p(z ) with a real argument z£l, which 
plays a predominant role in time-fractional diffusion, can also be efficiently computed with the Laplace transform 
and suitable quadrature rules (27) . 

The computation of the Wright function W P:P (z) is even more delicate. In theory, like before, it can be 
computed using power series for small values of the argument and a known asymptotic formula for large values, and 
for the intermediate case, values are obtained by using an integral representation [67]. The integral representation 
for the Wright function Wp, M (z) for intermediate values in the case of interest for the fundamental solution in 
one-dimension (where p = —a/2 < 0 , 0 </r = l + p<l and 2 = —x, x > 0 ) is given by 

POO 

w pA~ x )=* K(x,p,p,r)dr 
Jo 

where the kernel K{x, p, p, r ) is given by 

K(x, p, p,r) = r ~i 1 e ~ r + XCOB ( 7V P') r P sin(x sin(np)r~ p + tv p). 

This is a singular kernel with a leading order r~ p = r~ 1_p , with successive singular kernels of the form 
r -i- 3 p e -j- c , U p 0n expanding the terms. Hence, a direct treatment via numerical quadrature is inefficient, 
efficient approach is to use the change of variable s = r~ p , i.e., r — s _1 ^ p , and the transformed kernel is 

K(x, p, p, s) = (—p)~ 1 s < '~ p ^ (~p+ 1 'i~ 1 e - a p +xcoe(-7T p )s ^( 3 . s in( 7 rp)s + np). 

The fundamental solution of the one-dimensional time-fractional diffusion equation is expressed in terms of a 
Wright function W p , p (—x) with the choice p = —a/2 and p = 1 + p, cf. (12.211 . In this case the resulting kernel K 
simplifies to 

_ \ 

K(x, p , 1 + p, s) = {—p)~ 1 e~ s P + xcoa t*p') s sin(a;sin( 7 rp)s + (1 + p)io). 

This kernel is free from the grave singularity, and thus the quadrature method is quite effective. In general, the 
integral can be computed efficiently via the Gauss-Jacobi quadrature, with the weight function s ( '~ p ' > ( - ^+ 1 ) -1 . 

We note that an algorithm for the Wright function W p , p ,(z) over the whole complex plane C with rigorous error 
analysis is still missing. The endeavor in this direction would almost certainly involve dividing the complex 
domain C into different regions, and using different approximations on each region separately. 

A.2. Time fractional diffusion. We describe a finite difference method for the initial boundary value problem 
for the one-dimensional time-fractional diffusion equation 

dtU — u xx + qu = f(x, t) (x, t) £ fix ( 0 , T ), 

with the initial condition u(x, 0 ) = v and boundary conditions 

u(0,t) = g(t) and u(l,t) = h(t), t > 0 . 

There are many efficient numerical schemes for discretizing the problem. The discretization in space can be 
achieved by the standard central difference scheme, Galerkin finite element method m or spectral method, and 
the discretization in time can be achieved with the LI approximation |96ll63] and convolution quadrature [48]. We 
shall adopt the LI approximation in time and the central difference in space. Specifically, we divide the interval 
[0, T] into uniform subintervals, with nodes tk = hr, k = 0, ..., K, and a time step size t = T/K. Similarly, 
we partition the spatial domain fl into uniform subintervals, with grid points xt = ih, i = 0 ,..., N , and mesh 


r -i-2p, 

A more 
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size h = 1/N. Then the L 1-approximation of the Djrbashian-Caputo fractional derivative dtu(x,tk) developed 
in [96] [63] is given by: 


(A.l) 


d?u(x,t k ) 


— a 
T 


k -1 

bou(x , tk) ~ bk-iu(x, to) + ^2(bj - bj-i)u(x, tk-j) 
3=1 


where the weights bj are given by 

bj = (U + I) 1 ’ 0 - j 1 "“)/r (2 -a), j = 0,1 ,..., I< - 1 . 

If the solution u(x,t) is C 2 continuous in time, the local truncation error of the LI approximation is bounded 
by cr 2 ~ a for some c depending only on u [63l equation (3.3)]. In general, one can show that the scheme is only 
first-order accurate. Next with the central difference scheme in space and the notation Uj w u(xi,tk), we arrive 
at the following fully discrete scheme 


boUi - b k -!U° + y - bj-Ou'l J 

3=1 


. u i- 1 — 2U i + Ui + 1 k fk ■ i AT 1 

+ -^- + **=&, t = l,...,N-l, 


with Qi — g(xi) and f/ = f(xi, tk)- We note that at each time step, one needs to solve a tridiagonal linear system. 
However, the right hand side at the current step involves all previous steps, which can be quite expensive for a 
small step size, and this will most likely be required due to the first order in time convergence. This history piece 
represents one of the main computational challenges for time fractional differential equations. There are high-order 
schemes, e.g., convolution quadrature generated by the second-order backward difference formula [48]. Further, 
we note that the finite difference scheme in space can be replaced with the Galerkin finite element method, which 
is especially suitable for high dimensional problems on a general domain and elliptic operator involving variable 
coefficients m- 


A.3. Space fractional diffusion. Now we describe a fully discrete scheme based on the backward Euler method 
in time and a Galerkin finite element method in space for the space fractional diffusion problem on the unit 
interval fi = (0, 1) 

Ut — oD x u + qu = / infl x (0,7], 
with the initial condition u = v and the Dirichlet boundary condition 

u(0, t) — g(t) and u(l, t) = h(t), t > 0. 

The Galerkin finite element method relies on the variational formulation for the fractional elliptic problem 

qDx u + qu — f in fi, 

with a homogeneous Dirichlet boundary condition u{ 0) = u(l) = 0, recently developed in [45]. The variational 
formulation of the problem is given by: find u £ U = H^ 2 (Q.) such that 

~(oDx 12 u, xDi /2 v) + (qu, v) = (/, v) Vu £ V, 

where q D^v and are the left-sided and right-sided Riemann-Liouville derivative of order 7 £ (0,1) defined 

by (with c 7 = 1/T(1 — 7)) 

oDxv(x) = [ (x — s)~ 1 v(s)ds and ^Djv(x) — — C77— [ (x — s)~ 1 v(s)ds, 

dx J q dx J x 

respectively, and the test space V is given by 

V = jv £ H p/2 (fl) : u(l) = 0, (:e 1_/3 ,u) = oj . 

Now for the finite element discretization, we first divide the unit interval fi into a uniform mesh, with the grid 
points Xi = ih, i = 0,..., N and mesh size h = 1/N. Then for Uh C U we take the continuous piecewise linear 
finite element space, and for Vh C V we construct it from Uh- Specifically, with the finite element basis 4>i(x), 
i = 1,.. .,N — 1, we take 4>i(x) = (pi(x) — 7i(l — x) £ Vh, where the constant 7; is determined by the integral 
condition ( x4>i) = 0, i.e., 7j = h 2 ~^((i — l) 3- ' 9 + (i + 1) 3-/3 — 2i 3-/3 ). The computation of the leading term 
in the stiffness matrix and mass matrix can be carried analytically, and the part involving the potential q can be 
computed efficiently using quadrature rules; see [46] for details. 

Now for the time-dependent problem, like before, we divide the time interval [0, T\ into uniform subintervals, 
with tk = fer, k = 0,..., K, and the time step size r = T/K. Then with the backward Euler method in time, and 







INVERSE PROBLEMS FOR FDES 


27 


the finite element method in space, the approximate solution u £ at time tk can be split into u 1 ^ = + s k with 

the particular solution s k = g(tk) + ( h(tk ) — g(tk))x with the homogeneous solution £ Uh satisfying 

- (o ^Pi /2 v h ) + {quh,v h ) 

= ( f k , Vh) + r _1 (u^ _1 - - (gs\ %) Vvh £ 14. 

We note the resulting linear system is of lower Hessenberg form, due to the nonlocality of the fractional derivative 
operator. However, the coefficient matrix does not change during the time stepping procedure, and thus an LU 
factorization might be applied to speedup the computation. 
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